#include "mloop.h"
#include "ramachandran.h"

void init_Ramachandran_plot2 (char phi_psi[36][37], float phi_psi1[36][36][2])
{
	int i, j, k, n;
	char tmp;

	for (i = 0; i < 36; i++)
		for (j = 0; j < 36; j++)
		{
			tmp = phi_psi[i][j];
			phi_psi1[i][j][0] = phi_psi1[i][j][1] = 9999.0;
			if (tmp == 'x' || tmp == 'y' || tmp == 'z' || tmp == 'o')
			{
				for (n = 1; n <= 18; n++)
				{
					int p1, p2;

					p1 = p2 = 0;
					k = j + n;
					k = k % 36;
					tmp = phi_psi[i][k];
					if (tmp == 'x' || tmp == 'y' || tmp == 'z' || tmp == 'o')
						;
					else
					{
						p1 = 1;
						phi_psi1[i][j][0] = (k - 18) * 10.0 + 0.1;
					}

					k = j - n;
					k = (k + 36) % 36;
					tmp = phi_psi[i][k];
					if (tmp == 'x' || tmp == 'y' || tmp == 'z' || tmp == 'o')
						;
					else
					{
						p2 = 1;
						phi_psi1[i][j][1] = (k - 18)* 10.0 + 9.9;
					}
					if (p1 == 1 || p2 == 1)
						break;
				}
			}
			else
				phi_psi1[i][j][0] = phi_psi1[i][j][1] = -9999.0;
		}
}

void init_Ramachandran_plot3 (char phi_psi[36][37], float phi_psi2[36][36][2])
{
	int i, j, k, n;
	char tmp;

	for (i = 0; i < 36; i++)
		for (j = 0; j < 36; j++)
		{
			tmp = phi_psi[i][j];
			phi_psi2[i][j][0] = phi_psi2[i][j][1] = 9999.0;
			if (tmp == 'x' || tmp == 'y' || tmp == 'z' || tmp == 'o')
			{
				for (n = 1; n <= 18; n++)
				{
					int p1, p2;

					p1 = p2 = 0;
					k = i + n;
					k = k % 36;
					tmp = phi_psi[k][j];
					if (tmp == 'x' || tmp == 'y' || tmp == 'z' || tmp == 'o')
						;
					else
					{
						p1 = 1;
						phi_psi2[i][j][1] = (17 - k) * 10.0 + 9.9;
					}

					k = i - n;
					k = (k + 36) % 36;
					tmp = phi_psi[k][j];
					if (tmp == 'x' || tmp == 'y' || tmp == 'z' || tmp == 'o')
						;
					else
					{
						p2 = 1;
						phi_psi2[i][j][0] = (17 - k) * 10.0 + 0.1;
					}
					if (p1 == 1 || p2 == 1)
						break;
				}
			}
			else
				phi_psi2[i][j][0] = phi_psi2[i][j][1] = -9999.0;
		}
}

void init_Ramachandran_plot ()
{
	init_Ramachandran_plot2 (phi_psi, phi_psi1);
	init_Ramachandran_plot3 (phi_psi, phi_psi2);
	init_Ramachandran_plot2 (phi_psi_pro, phi_psi_pro1);
	init_Ramachandran_plot3 (phi_psi_pro, phi_psi_pro2);
	init_Ramachandran_plot2 (phi_psi_gly, phi_psi_gly1);
	init_Ramachandran_plot3 (phi_psi_gly, phi_psi_gly2);
}

/* normal angle from range (-360, 360) to rangle (-180, 180) */
float normal_angle(float angle)
{
	if (angle > 180.0)
		angle = angle - 360.0;
	else if (angle < -180.0)
		angle = angle + 360.0;

	return angle;
}

float allow_angle2(float phi, float psi, float Ramachandran[36][36][2], float angle)
{
	int i, j, k, n;
	float theta, angle0, angle1, tmp;

	i = (phi + 180.0) / 10.0;
	j = (psi + 180.0) / 10.0;
	j = 35 - j;
	tmp = Ramachandran[j][i][0];
	if (tmp < -1000.0)
		theta = 0.0;
	else if (Ramachandran[j][i][0] > 1000.0 && Ramachandran[j][i][1] > 1000.0)
		theta = 0.0;
	else
	{
		angle0 = normal_angle (Ramachandran[j][i][0] - angle);
		angle1 = normal_angle (Ramachandran[j][i][1] - angle);
		if (Ramachandran[j][i][0] > 1000.0) 
			theta = angle1;
		else if (Ramachandran[j][i][1] > 1000.0)
			theta = angle0;
		else if (fabs (angle0) > fabs (angle1))
			theta = angle1;
		else
			theta = angle0;
	}

	return theta;
}

float allow_angle(float phi, float psi, int type, int flag, float angle)
{
	float theta;

	if (flag == 0)
	{
		phi += angle;
		phi = normal_angle (phi);
	}
	else
	{
		psi += angle;
		psi = normal_angle (psi);
	}

	if (type == 6 && flag == 0)
		theta = allow_angle2 (phi, psi, phi_psi_gly1, phi);
	else if (type == 6 && flag == 1)
		theta = allow_angle2 (phi, psi, phi_psi_gly2, psi);
	else if (type == 12 && flag == 0)
		theta = allow_angle2 (phi, psi, phi_psi_pro1, phi);
	else if (type == 12 && flag == 1)
		theta = allow_angle2 (phi, psi, phi_psi_pro2, psi);
	else if (flag == 0)
		theta = allow_angle2 (phi, psi, phi_psi1, phi);
	else
		theta = allow_angle2 (phi, psi, phi_psi2, psi);

	return theta;
}

/* phi, psi (-180, 180) */
int judge_phi_psi(float phi, float psi, int type)
{
	int i, j;
	char flag;

	i = (phi + 180.0) / 10.0;
	j = (psi + 180.0) / 10.0;

	/* gly */
	if (type == 6)
		flag = phi_psi_gly[35 - j][i];
	/* pro */
	else if (type == 12)
		flag = phi_psi_pro[35- j][i];
	else
		flag = phi_psi[35 - j][i];

	if (flag == 'x' || flag == 'y' || flag == 'z' || flag == 'o')
		return 0;
	else
		return 1;
}

void rand_phi_psi(float *phi, float *psi, int type)
{
	float a, b;

	for (;;)
	{
		a = 360.0 * drand48();
		b = 360.0 * drand48();
		a = a - 180.0;
		b = b - 180.0;
		if (judge_phi_psi (a, b, type))
			break;
	}
	*phi = a;
	*psi = b;
}

void rand_phi_psi2(float *phi, float *psi, float phi_ref, float psi_ref, int type, float range)
{
	float a, b;

	for (;;)
	{
		a = 2.0 * range * drand48();
		a = phi_ref + a - range;
		a = normal_angle (a);

		b = 2.0 * range * drand48();
		b = psi_ref + b - range;
		b = normal_angle (b);

		if (judge_phi_psi (a, b, type))
			break;
	}
	*phi = a;
	*psi = b;
}
